/*
 * Cell.cpp
 *
 *  Created on: Feb 28, 2011
 *      Author: wjfortin
 */

#include "Cell.h"
#include <string.h>
#include <sstream>

int Cell::cellCounter = 0;

/**
 * Cell constructor. 
 *
 * The Cell contains all information about its kinetic equations (in the form of a Directed Graph) describing the concentration of the Molecules it contains.
 * The Cell constructor also creates the directory structure to hold associated output files. 
 *
 */
Cell::Cell() {

	//give each cell a unique ID
	cellID = cellCounter++;

	numGenerations = 0;
	
	//Directed graph which describes the interactions between molecules.
	equations = new ListDigraph();

	moleculeCount = 0;
	
	numProtein = 0;
	numMRNA = 0;
	numDNA = 0;

	//Map Molecule objects onto vertices
	vertexes = new ListDigraph::NodeMap<Molecule*>(*equations);

	//Map Interaction objects onto edges
	edges = new ListDigraph::ArcMap<Interaction*>(*equations);

	//the nullNode is the target of all degredation interactions
	nullNode = addMRNA(0.0);
	
	//the concentration of degraded molecules need not be calculated
	(*vertexes)[nullNode]->setCalculate(false);
	
	//initialize the starting molecules in the cell
	ListDigraph::Node A = addProtein(DEFAULT_PROTEIN_CONC);
	ListDigraph::Node B = addMRNA(DEFAULT_MRNA_CONC);
	ListDigraph::Node C = addDNA(DEFAULT_DNA_CONC);
	
	//get random rates to assign to the kinetic interactions
	double kTXN = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	double kTSLN = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	double kDEGp = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	double kDEGmRNA = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	
	//create interactions between the initial molecules with the above rates
	addInteraction(B,A, kTSLN, 0.0, 1.0, "tsln");
	addInteraction(C,B, kTXN, 0.0, 1.0, "txn");
	addInteraction(A,nullNode,kDEGp, 1.0, 0,"deg");
	addInteraction(B,nullNode,kDEGmRNA, 1.0, 0,"deg");

	//create directory structure for output
	char basedir[80];
	char celldir[80];
	char mutationFile[80];

	//create base folder name with the current process id
	sprintf(basedir, "./output/%d", getpid());
	//subfolder name for this cell
	sprintf(celldir, "./output/%d/cell%d", getpid(), cellID);
	//filename for mutation log
	sprintf(mutationFile, "./output/%d/cell%d/mutations.log", getpid(), cellID);

	//create the directories
	mkdir("./output", S_IRWXU | S_IRWXG | S_IRWXO);
	mkdir(basedir, S_IRWXU | S_IRWXG | S_IRWXO);
	mkdir(celldir, S_IRWXU | S_IRWXG | S_IRWXO);



	//open mutation log file
	mutationLog.open(mutationFile);

	//error check for opening mutationLog file
	if(!mutationLog)
		perror("mutationLog");


	TRACE(2,"Cell " << cellID << " created");
}
/**
 * Cell destructor.
 * 
 * Close output file, delete molecule and interaction objects, and delete other containers.
 */
Cell::~Cell() {

	// flush pending output and close the file cleanly
	mutationLog.flush();
	mutationLog.close();

	TRACE(3, "Mutation log flushed.");
	
	// delete the Molecule objects from the graph
	for (ListDigraph::NodeIt node(*equations); node != INVALID; ++node){
		delete[] (*vertexes)[node]->rkData;
		delete (*vertexes)[node];
	}
	
	// delete the Interaction objects from the graph
	for (ListDigraph::ArcIt edge(*equations); edge != INVALID; ++edge)
		delete ((*edges)[edge]);

	// clear the proteinList container
	proteinList.clear();

	// delete the graph objects
	delete vertexes;
	delete edges;
	delete equations;
}

/**
 * Increment the current generation number, and mutate the Cell. 
 *
 * @return The return value of the mutate() method
 */
int Cell::nextGeneration()
{
	numGenerations++;
	return mutate();
}

/**
 * Get a random Node which contains a Protein. 
 *
 * The proteinList vector holds references to proteins which are added when a new protein is created. A random integer is used to index into the vector. 
 *
 * @return A random Node object in the graph containing a Protein.
 */
ListDigraph::Node Cell::getRandomProtein()
{
	TRACE(8,"Getting random number...");

	//select a random index from the protein list
	int num = r.randInt(proteinList.size()-1);

	TRACE(8,"Got random number...");

	//return the node containing that protein
	return equations->nodeFromId(proteinList.at(num)->nodeID);

}

/**
 * Add a protein with the given initial concentration to the Cell.
 *
 * @param concentration Initial concentration for the new Protein
 *
 * @return A Node object in the graph containing a newly created Protein.
 */
ListDigraph::Node Cell::addProtein(double concentration)
{

	TRACE(4,"Adding Protein in Cell " << cellID << "...");
	
	//add a node to the graph
	ListDigraph::Node newNode = equations->addNode();
	//create a new Protein object in the node
	(*vertexes)[newNode] = new Protein();
	//with the nodeID, a Protein can find its location within the graph
	(*vertexes)[newNode]->nodeID = equations->id(newNode);
	//set concentration
	(*vertexes)[newNode]->getConcentration()->value = concentration;
	//set initial concentration
	(*vertexes)[newNode]->getInitialConcentration()->value = concentration;

	//add protein to the protein list
	proteinList.push_back((Protein*)(*vertexes)[newNode]);

	TRACE(4,"Protein in Cell " << cellID << " added.");
	//increment molecule counter
	moleculeCount++;
 	numProtein++;
 	
	return newNode;
}

/**
 * Add an existing protein object to the Cell.
 *
 * @param protein Pointer to a protein to add to the graph.
 *
 * @return A Node object in the graph containing the newly added Protein.
 */
ListDigraph::Node Cell::addProtein(Protein* protein)
{

	TRACE(4,"Adding Protein in Cell " << cellID << "...");
	//add a node to the graph
	ListDigraph::Node newNode = equations->addNode();
	//store the protein object in the newly created node
	(*vertexes)[newNode] = protein;
	//with the nodeID, a Protein can find its location within the graph
	(*vertexes)[newNode]->nodeID = equations->id(newNode);
	//set initial concentration
	(*vertexes)[newNode]->getInitialConcentration()->value = (*vertexes)[newNode]->getConcentration()->value;

	//add protein to the protein list
	proteinList.push_back((Protein*)(*vertexes)[newNode]);

	TRACE(4,"Protein in Cell " << cellID << " added.");
	//increment molecule counter
	moleculeCount++;
	numProtein++;

	return newNode;
}

/**
 * Add an MRNA with the given initial concentration to the Cell.
 *
 * @param concentration Initial concentration for the new MRNA
 *
 * @return A Node object in the graph containing a newly created MRNA.
 */
ListDigraph::Node Cell::addMRNA(double concentration)
{

	TRACE(4,"Adding MRNA in Cell " << cellID << "...");
	ListDigraph::Node newNode = equations->addNode();
	(*vertexes)[newNode] = new MRNA();
	(*vertexes)[newNode]->nodeID = equations->id(newNode);
	(*vertexes)[newNode]->getConcentration()->value = concentration;
	(*vertexes)[newNode]->getInitialConcentration()->value = concentration;

	TRACE(4,"MRNA in Cell " << cellID << " added.");
	moleculeCount++;
	numMRNA++;
	
	return newNode;
}

/**
 * Add a DNA with the given initial concentration to the Cell.
 *
 * @param concentration Initial concentration for the new DNA
 *
 * @return A Node object in the graph containing a newly created DNA.
 */
ListDigraph::Node Cell::addDNA(double concentration)
{

	TRACE(4,"Adding DNA in Cell " << cellID << "...");
	ListDigraph::Node newNode = equations->addNode();
	(*vertexes)[newNode] = new DNA();
	(*vertexes)[newNode]->nodeID = equations->id(newNode);
	(*vertexes)[newNode]->getConcentration()->value = concentration;
	(*vertexes)[newNode]->getInitialConcentration()->value = concentration;

	TRACE(4,"DNA in Cell " << cellID << " added.");
	moleculeCount++;
	numDNA++;
	
	return newNode;
}

/**
 * Add an interation between two nodes in the Cell.
 *
 * @param from The Node in the graph containing the substrate Molecule
 * @param to The Node in the graph containing the product Molecule
 * @param rate The kinetic rate at which the reaction takes place
 * @param substrateFactor A constant multiplied factor of the amount of substrate Molecule consumed
 * @param productFactor A constant multiplied factor of the amount of product Molecule produced
 * @param interactionName A description for the interaction (fwd, rev, deg, txn, tsln)
 *
 * @return An Arc object in the graph containing the newly created Interaction.
 */
ListDigraph::Arc Cell::addInteraction(ListDigraph::Node from, ListDigraph::Node to, double rate, double substrateFactor, double productFactor, const char* interactionName){

	TRACE(4,"Adding Interaction in Cell " << cellID << "...");
	ListDigraph::Arc newArc = equations->addArc(from,to);
	Quantity* interactionRate = new Quantity();
	interactionRate->value = rate;
	(*edges)[newArc] = new Interaction(interactionRate, substrateFactor, productFactor, interactionName);

	TRACE(4,"Interaction in Cell " << cellID << " added.");
	return newArc;
}

/**
 * Add an interation between two nodes in the Cell.
 *
 * @param from The Node in the graph containing the substrate Molecule
 * @param to The Node in the graph containing the product Molecule
 * @param rate A pointer to a Quantity object containing the rate at which the reaction takes place 
 * @param substrateFactor A constant multiplied factor of the amount of substrate Molecule consumed
 * @param productFactor A constant multiplied factor of the amount of product Molecule produced
 * @param interactionName A description for the interaction (fwd, rev, deg, txn, tsln)
 *
 * @return An Arc object in the graph containing the newly created Interaction.
 */
ListDigraph::Arc Cell::addInteraction(ListDigraph::Node from, ListDigraph::Node to, Quantity* rate, double substrateFactor, double productFactor, const char* interactionName){

	TRACE(4,"Adding Interaction in Cell " << cellID << "...");
	ListDigraph::Arc newArc = equations->addArc(from,to);
	Quantity* interactionRate = new Quantity();
	interactionRate = rate;
	(*edges)[newArc] = new Interaction(interactionRate, substrateFactor, productFactor, interactionName);

	TRACE(4,"Interaction in Cell " << cellID << " added.");
	return newArc;
}

/**
 * Select a random mutation to carry out on the cell during the current generation.
 *
 * The mutate method randomly chooses a mutation type. The effect of the mutation is randomly chosen, and applied to the Cell.  
 * Mutation types:
 *  Forward Rate Change - Modify the rate of a randomly chosen Interaction (Interaction::interactionName = 'fwd' || 'tsln' || 'txn')
 *  Reverse Rate Change - Modify the rate of a randomly chosen Interaction (Interaction::interactionName = 'rev')
 *  Degradation Rate Change - Modify the rate of a randomly chosen Interaction (Interaction::interactionName = 'deg')
 *  New Post Translational Modification - Randomly choose a Protein, and add Interactions for the production of a new PTM
 *  Histone Modification - Modify the productFactor of a randomly chosen Interaction (Interaction::interactionName = 'tsln')
 *  New Complex - Randomly choose two Proteins, and add Interactions for them to form a Complex together
 *  New Basic Protein - Create a new Basic Protein
 *  New Protein-Promoter Complex - Not Implemented
 *
 *  The Mutation is logged to the mutation.log file in the output folder for the Cell.
 *
 *  @return (Unused/Unimplemented) Return code describing success of Mutation
 */
int Cell::mutate()
{
	TRACE(3,"Mutating Cell " << getCellID() << "...");
	// add a new line to the mutation log
	mutationLog << "\nGeneration " << numGenerations << endl;
	
	// pick category of mutation from weighted probabilities: large, small, or null
	double category = r.rand(1.0);

	// pick subcategory from weighted probabilities
	double subcategory = r.rand(1.0);

	int result = 0;

	//TODO 6/4/11 remove magic numbers
     //small mutations
     // TODO refactor mutation type probabilities to Constants.h 
	if (category < .4) {
		if (subcategory < .20)
			result = fwdRateChange();
		else if (subcategory < .40)
			result = revRateChange();
		else if (subcategory < .60)
			result = degRateChange();
		else if (subcategory < .80)
			result = newPostTransMod();
		else
			result = histoneMod();
	}
	//large mutations
	else if (category < .7) {
		if (subcategory < .33)
			result = newComplex();
		else if (subcategory < .67)
			result = newBasic();
		else
			result = newPromoter();
	}
	//null mutation
	else {
		mutationLog << " Null mutation" << endl;
		result = 0;
	}

	TRACE(2,"Cell " << getCellID() << " mutated.");

	return result;
}

/**
 * Return the current CellID.
 *
 * Upon creation, each Cell receives a unique increasing ID.
 *
 * @return Unique identifier for the current Cell.
 */
int Cell::getCellID()
{
	return cellID;
}



/**
 * Output the data generated by Runge-Kutta for all molecules for the current generation
 *
 * Saves the Runge-Kutta data for the current generation in output files.
 * The current configuration of the Cell, showing all interactions, rates, and factors between Molecules is also logged.
 *
 * After Cell::rungeKuttaEvaluate() has been called, each Molecule holds an internal array (rkData) with the data points calculated by Runge-Kutta.
 * Each Molecule in the cell writes data to an output file for the current cell and generation.
 *
 * If the data was produced by Runge-Kutta using the RUNGE_KUTTA_PRECISE_STEP / RUNGE_KUTTA_PRECISE_LIMIT parameters, the following additional actions take place:
 *
 * The data is plotted to png files by Gnuplot. Each Molecule in the Cell will produce it's own graph, as well as a single graph with all Molecules plotted.
 *
 */
void Cell::outputRungeKuttaData(){

	TRACE(1, "Beginning display...");

	Molecule* currentNode;
	Interaction* currentEdge;

	ofstream output;

	char filename[181];

	//create the output file name for this cell
	// the file will be named based on the process id (to be unique) the cell number, and the runge-kutta precision
	// e.g.   1043.cell0.std     1043.cell1.std     1043.cell1.precise
	if(rkDataStep == RUNGE_KUTTA_STD_STEP && rkDataLimit == RUNGE_KUTTA_STD_LIMIT)
		snprintf(filename, 179, "./output/%d/cell%d/gen%d.std", getpid(), cellID, numGenerations);

	else if(rkDataStep == RUNGE_KUTTA_PRECISE_STEP && rkDataLimit == RUNGE_KUTTA_PRECISE_LIMIT)
		snprintf(filename, 179, "./output/%d/cell%d/gen%d.precise", getpid(), cellID, numGenerations);
	else
		snprintf(filename, 179, "./output/%d.cell%d.gen%d", getpid(), cellID, numGenerations);

	//open the output file
	output.open(filename);
	//cout << filename << endl;
	if(!output)
		perror("output");
	//cout << "error";

	//output a rough configuration for the cell (Will be improved / cleaned up)
	output << "Cell " << cellID << " configuration:\n";

	// display all interactions between the molecules and their k values
	for (ListDigraph::ArcIt a(*equations); a != INVALID; ++a){
		currentEdge = (*edges)[a];

		Molecule* source = (*vertexes)[equations->source(a)];
		Molecule* target = (*vertexes)[equations->target(a)];

		output << "[";

		if(source->calculate == true)
			output << source->name << source->getMoleculeID();
		else
			output << "null";

		output << "]--" << currentEdge->getSubstrateFactor() << "--(" << currentEdge->getRate() << " )--" << currentEdge->getProductFactor() << "-->[";

		if(target->calculate == true)
			output << target->name << target->getMoleculeID();
		else
			output << "null";

		output << "]\n";

	}


	// rkDataLimit and rkDataStep describe the limit and stepsize of the last runge-kutta evaluation
	// dividing the limit by the stepsize gives us the number of points in each array
	int arraySize = (int) ceil(rkDataLimit/rkDataStep);
        
        FILE* completePlot;
	
	// if this is a precise run, open a graph file
	// completePlot will plot all molecules to a single graph
	if(rkDataStep == RUNGE_KUTTA_PRECISE_STEP && rkDataLimit == RUNGE_KUTTA_PRECISE_LIMIT){
	        
		completePlot = popen("gnuplot > /dev/null 2>&1","w");

	        fprintf(completePlot,"set term png\n");
	        fflush(completePlot);

                fprintf(completePlot,"set xlabel \"time\"\n");
                fflush(completePlot);

                fprintf(completePlot,"set ylabel \"conc\"\n");
                fflush(completePlot);

                fprintf(completePlot,"set format x \"%c03.2f\"\n",'%');
                fflush(completePlot);

	        //set the filename of the graph
	        fprintf(completePlot,"set output \"%s.complete.png\"\n",filename);
		fflush(completePlot);
		
		fprintf(completePlot,"plot ");
		fflush(completePlot);

	  }
	// iterate through all of the molecules
	for(ListDigraph::NodeIt n(*equations); n!=INVALID;++n){
		currentNode = (*vertexes)[n];

		if(!currentNode->calculate)
			continue;

		//output the runge-kutta values for that molecule to the cells output file
		for(int i = 0; i < arraySize; i++){
			TRACE(7,currentNode->getMoleculeID() << " - " << i << " " << currentNode->rkData[i]);
			output << currentNode->getMoleculeID() << " - " << i*rkDataStep << " " << currentNode->rkData[i] << endl;
		}
	
		//only generate graphs for precise RK output
		if(rkDataStep != RUNGE_KUTTA_PRECISE_STEP || rkDataLimit != RUNGE_KUTTA_PRECISE_LIMIT)
		  continue;

		//do not generate graphs for DNA molecules
		if(!strcmp(currentNode->name.c_str(),"DNA"))
		  continue;
	
        	FILE *gnuplot;
	    
	        //open a gnuplot process and get a pipe to its stdin
                gnuplot = popen("gnuplot","w");
		
	        //set png format
	        fprintf(gnuplot,"set term png\n");
	        fflush(gnuplot);

                fprintf(gnuplot,"set xlabel \"time\"\n");
                fflush(gnuplot);

                fprintf(gnuplot,"set ylabel \"conc\"\n");
                fflush(gnuplot);

                fprintf(gnuplot,"set format x \"%c03.2f\"\n",'%');
                fflush(gnuplot);

	        //set the filename of the graph
	        fprintf(gnuplot,"set output \"%s.mol%d.png\"\n",filename, currentNode->getMoleculeID() );
		fflush(gnuplot);

                //plot to individual file
                fprintf(gnuplot,"plot \"%s\" using 3:($1==%d ? $4 : 1/0) t \"Molecule %d\"\n",filename, currentNode->getMoleculeID(), currentNode->getMoleculeID());
                fflush(gnuplot);

		//end the gnuplot process
                fprintf(gnuplot,"exit\n");
                fflush(gnuplot);

		//close the process pipe
                pclose(gnuplot);


		// plot to the complete graph file
                fprintf(completePlot,"\"%s\" using 3:($1==%d ? $4 : 1/0) t \"Molecule %d\" pt 1 , ",filename, currentNode->getMoleculeID(), currentNode->getMoleculeID());
		fflush(completePlot);

	
	}

	if(rkDataStep == RUNGE_KUTTA_PRECISE_STEP && rkDataLimit == RUNGE_KUTTA_PRECISE_LIMIT)
	{	
		// terminate the plot command cleanly
                fprintf(completePlot,"\"%s\" using 3:($1==-1 ? 1/0 : 1/0) t \"\" pt 1 \n ",filename);
		fflush(completePlot);

		//end the gnuplot process
		fprintf(completePlot,"exit\n");
		fflush(completePlot);

		//close the process pipe
		pclose(completePlot);
	}
	//close the output file
	output.close();

	TRACE(1, "Displayed.");

}

/**
 * Assign a score to the Cell based on Runge-Kutta data calculated for Molecules. 
 *
 * The scoring method checks successive data points, and counts the number of times the slope changes between positive or negative.
 * 
 * The Cell's score is then taken to be equal to the Molecule with the highest score (most number of changes).
 *
 * @return The score of the Runge-Kutta data of the current Cell/Generation as determined by the evaluation function.
 */
int Cell::getScore(int scoringType)
{

	int score = -1;

	// use the chosen scoring function
        switch(scoringType){
	
	//default behavior
	default:
	//scoring function 1
	case 1:
		Molecule* currentNode;


		int prevDir = 0;
		int currentDir = 0;
		int numChanges = 0;
		int maxChange = 0;

     		// get the number of data points held in rkData array
     		// the time limit divided by the step size is equivalent to the total number of steps, or data points
		int arraySize = (int) ceil(rkDataLimit/rkDataStep);

		// for each molecule in the cell
		for(ListDigraph::NodeIt n(*equations); n!=INVALID;++n){
			currentNode = (*vertexes)[n];

			prevDir = 0;
			currentDir = 0;
			numChanges = 0;

			// count the number of sign changes in the rkData
			for(int j = 1; j< arraySize; j++){

				//value has stayed constant, move to next point
				if(currentNode->rkData[j] - currentNode->rkData[j-1] <= .01 && currentNode->rkData[j] - currentNode->rkData[j-1] >= -.01)
					continue;

				//value has decreased
				if(currentNode->rkData[j] - currentNode->rkData[j-1] < -.01)
					currentDir = -1;

				//value has increased
				if(currentNode->rkData[j] - currentNode->rkData[j-1] > .01) 
					currentDir = 1;


				// now increasing, was decreasing
				if( currentDir == 1 && prevDir == -1)
					numChanges ++;

				// now decreasing, was increasing
				if( currentDir == -1 && prevDir == 1)
					numChanges ++;

				//store new direction
				prevDir = currentDir;

			}

			//keep max
			if(numChanges > maxChange)
				maxChange = numChanges;

		}
		//assign the score to return
		score = maxChange;


	}// end switch statement

	//send data to output files for this generation
	outputRungeKuttaData();

	//should not ever reach here
	return score;
}

/**
 * Calculate concentration values for molecules in the cell over time, using the 4th order Runge-Kutta method.
 * 
 * The Directed Graph (ListDigraph equations) is representative of the differential equations which describe the concentrations of the Molecules. This implementation of
 * Runge-Kutta operates directly on the DiGraph to calculate successive data points (x,y) = (t*rkTimestep, rkData[t]) while t*rkTimestep < rkTimeLimit
 */
void Cell::rungeKuttaEvaluate(double rkTimestep, double rkTimeLimit) {
  
     // save the step and limit size for later
     // when using rkData, knowing the step and limit is important to distinguish between standard precision / higher precision data, 
     //as well as knowing how many elements are in the array
	rkDataStep = rkTimestep;
	rkDataLimit = rkTimeLimit;
	
	//cleanRkData allocates enough space for the rkData array in each molecule to hold the data that Runge-Kutta is about to produce
	cleanRkData(rkTimestep, rkTimeLimit);
	
	int t = 0;
	//initial t0 for numerical integration
	double currentTime = 0;

     // temporary pointers 
	Molecule* currentNode;
	Molecule* oppositeNode;
	Interaction* currentEdge;

     // loop through each molecule, and reset each molecule to begin evaluation at its given initial concentration
	for(ListDigraph::NodeIt n(*equations); n != INVALID; ++n){
		currentNode = (*vertexes)[n];

		// skip the null node
		if(!currentNode->calculate)
			continue;

		TRACE(6, "rkdata["<<(int)ceil(rkTimeLimit/rkTimestep)<<"]");
		//update the concentration to initial concentration
		currentNode->getConcentration()->value = currentNode->getInitialConcentration()->value;
         //initialize the first rkData point to the molecules initial concentration
		currentNode->rkData[0] = currentNode->getConcentration()->value;

		TRACE(6, currentNode->getMoleculeID() << " rkData[0]" << currentNode->getConcentration()->value);
	}

     // while Runge-Kutta still has points to calculate
	while(t < ceil(rkTimeLimit/rkTimestep)){

		TRACE(6,"Time = " << currentTime << " (limit " << rkTimeLimit << ")");


		//iterate through the 4 intermediate k values
		for(int k = 0; k<4; k++){

			TRACE(6,"k = " << k);

			//for each node in the graph
			for(ListDigraph::NodeIt n(*equations); n!=INVALID;++n){
				currentNode = (*vertexes)[n];

				// skip the null node
				if(!currentNode->calculate)
					continue;

				//for each outgoing edge from that node
				for(ListDigraph::OutArcIt a(*equations,n); a!=INVALID;++a){
					currentEdge = (*edges)[a];

					TRACE(6,"Cell " << cellID << " Molecule " << currentNode->getMoleculeID() << " -> molecule "<< (*vertexes)[(*equations).oppositeNode(n,a)]->getMoleculeID())
					TRACE(6,"Cell " << cellID << " Molecule " << currentNode->getMoleculeID() << " old rk[" << k << "]="<<currentNode->getConcentration()->rkValues[k])

					//update the k value for runge-kutta calculation
					switch (k){
					case 0:
						currentNode->getConcentration()->rkValues[0] -= currentEdge->getSubstrateFactor() * currentEdge->getRate(0) * currentNode->getConcentration()->value;
						break;
					case 1:
						currentNode->getConcentration()->rkValues[1] -= currentEdge->getSubstrateFactor() * (currentEdge->getRate(1) * (currentNode->getConcentration()->value + (rkTimestep/2 * currentNode->getConcentration()->rkValues[0])));
						break;
					case 2:
						currentNode->getConcentration()->rkValues[2] -= currentEdge->getSubstrateFactor() * (currentEdge->getRate(2) * (currentNode->getConcentration()->value + (rkTimestep/2 * currentNode->getConcentration()->rkValues[1])));
						break;
					case 3:
						currentNode->getConcentration()->rkValues[3] -= currentEdge->getSubstrateFactor() * (currentEdge->getRate(3) * (currentNode->getConcentration()->value + (rkTimestep * currentNode->getConcentration()->rkValues[2])));
						break;

					}
					TRACE(6,"Cell " << cellID << " Molecule " << currentNode->getMoleculeID() << " new rk[" << k << "]="<<currentNode->getConcentration()->rkValues[k])

				}
				//for each incoming edge to that node
				for(ListDigraph::InArcIt a(*equations,n); a!= INVALID;++a){
					currentEdge = (*edges)[a];
					oppositeNode = (*vertexes)[(*equations).oppositeNode(n,a)];

					TRACE(6,"Cell " << cellID << " Molecule " << currentNode->getMoleculeID() << " <- molecule "<< oppositeNode->getMoleculeID())

					TRACE(6,"Cell " << cellID << " Molecule " << currentNode->getMoleculeID() << " old rk[" << k << "]="<<currentNode->getConcentration()->rkValues[k])

					//update the k value for runge kutta calculation
					switch (k){
					case 0:
						currentNode->getConcentration()->rkValues[0] += currentEdge->getRate(0) * currentEdge->getProductFactor() * oppositeNode->getConcentration()->value;
						break;
					case 1:
						currentNode->getConcentration()->rkValues[1] += currentEdge->getRate(1) * currentEdge->getProductFactor() * (oppositeNode->getConcentration()->value + (rkTimestep/2 * oppositeNode->getConcentration()->rkValues[0]));
						break;
					case 2:
						currentNode->getConcentration()->rkValues[2] += currentEdge->getRate(2) * currentEdge->getProductFactor() * (oppositeNode->getConcentration()->value + (rkTimestep/2 * oppositeNode->getConcentration()->rkValues[1]));
						break;
					case 3:
						currentNode->getConcentration()->rkValues[3] += currentEdge->getRate(3) * currentEdge->getProductFactor() * (oppositeNode->getConcentration()->value + (rkTimestep * oppositeNode->getConcentration()->rkValues[2]));
						break;
					}
					TRACE(6,"Cell " << cellID << " Molecule " << currentNode->getMoleculeID() << " new rk[" << k << "]="<<currentNode->getConcentration()->rkValues[k])

				}
				TRACE(6,"Molecule " << currentNode->getMoleculeID() << " rkVal[" << k << "] = " << currentNode->getConcentration()->rkValues[k]);

			}
		}

		//advance the timestep
		currentTime += rkTimestep;
		t += 1;

		//update the concentration values for the molecules in the cell
		for(ListDigraph::NodeIt n(*equations); n!=INVALID;++n){
			currentNode = (*vertexes)[n];

			//skip the null node
			if(!currentNode->calculate)
				continue;
              // calculate the new concentration, based on the previous concentration and the four calculated rkValues
			double newConcentration = ( currentNode->rkData[t - 1] + ((rkTimestep/6))*( currentNode->getConcentration()->rkValues[0] + 2 * currentNode->getConcentration()->rkValues[1] + 2 * currentNode->getConcentration()->rkValues[2] + currentNode->getConcentration()->rkValues[3]));
			
			// a negative concentration is not possible in this context
			if (newConcentration < 0)
				newConcentration = 0;
				
			// reset the rkValues to 0 for next iteration
			for (int i = 0; i < 4; i++)
				currentNode->getConcentration()->rkValues[i] = 0;

			TRACE(6,currentNode->getMoleculeID() << " rkData[" << t << "] = " << newConcentration);

              // save the newly calculated concentration in the molecule's rkData array
			currentNode->rkData[t] = newConcentration;
			// update the current concentration of the molecule
			currentNode->getConcentration()->value = newConcentration;

		}
	}
}

/**
 * Set up the rkData array in each molecule.
 *
 * The rkData array holds the concentration data calculated by Runge-Kutta. Plotted as (x,y) = (t*step,rkData[t])
 * rkData needs to be large enough to hold (limit / step) data points. To ensure this rkData is deleted and reallocated before each call to Runge-Kutta.
 */
void Cell::cleanRkData(double step, double limit){


	Molecule* currentNode;
	//loop through all molecules in the cell
	for(ListDigraph::NodeIt n(*equations); n!=INVALID;++n)
	{
		currentNode = (*vertexes)[n];
          
         // delete the old rkData array
		delete[] currentNode->rkData;
         
         // create rkData of the correct size
		currentNode->rkData = new double[ (int) ceil(limit/step) + 1];
		//initialize the values in the array to 0
		for (int i = 0; i < (int) ceil(limit/step) + 1; i++)
		{
			currentNode->rkData[i] = 0;
		}

	}
}


/**
 * Fill in
 */

bool Cell::isComplexed(Protein* p1, Protein* p2){
	
	TRACE(5, "Checking for pre-existing complex between Proteins " << p1->getMoleculeID() << " and " << p2->getMoleculeID() << "...");
	Molecule* oppositeNode1;
	Molecule* oppositeNode2;

	ListDigraph::Node n1 = equations->nodeFromId(p1->nodeID);
	ListDigraph::Node n2 = equations->nodeFromId(p2->nodeID);

	//	ListDigraph::NodeIt n1(*equations, a1);
	//	ListDigraph::NodeIt n2(*equations, a2);

	// for each protein on an outgoing edge from the first protein's node...
	for (ListDigraph::OutArcIt a1(*equations, n1); a1 != INVALID; ++a1){
		oppositeNode1 = (*vertexes)[(*equations).oppositeNode(n1,a1)];

		// check if the second protein has an outgoing edge to the same protein
		for (ListDigraph::OutArcIt a2(*equations, n2); a2 != INVALID; ++a2){

			oppositeNode2 = (*vertexes)[(*equations).oppositeNode(n2,a2)];

			TRACE(6,oppositeNode1->getMoleculeID());

			if (oppositeNode1->getMoleculeID() == oppositeNode2->getMoleculeID() && typeid(*oppositeNode1) == typeid(Protein)){

				TRACE(5, "Pre-existing complex found: " << oppositeNode1->name << oppositeNode1->getMoleculeID());
				TRACE(5, "Protein 1: " << p1->name << p1->getMoleculeID() << typeid(*p1).name() << " and Protein 2: " << p2->name << p2->getMoleculeID() << typeid(*p2).name());
				return true;
			}
		}
	}

	TRACE(5, "Pre-existing complex not found");
	return false;
}

/**
 * Fill in
 */
int Cell::newComplex()
{
	/*
	//placeholder message
	mutationLog << " New Complex (not implemented)" << endl;
	return -1;
	 */

	mutationLog << " New Complex" << endl;

	TRACE(3,"Creating Protein-Protein Complex in Cell " << cellID << "...");

	ListDigraph::Node n1 = getRandomProtein();
	ListDigraph::Node n2 = getRandomProtein();
	Protein* p1 = (Protein*) (*vertexes)[n1];
	Protein* p2 = (Protein*) (*vertexes)[n2];

	//If there is no limit on the number of tries, the loop below could go forever.
	//The limit is based on the number of proteins; currently this limit is twice
	//the number of proteins.
	int tries = proteinList.size() * 2;

	//makes sure at least one try is made
	tries = tries <= 1 ? tries : 1;

	while ((p1->getMoleculeID() == p2->getMoleculeID() || isComplexed(p1, p2)) && tries > 0) {

		n1 = getRandomProtein();
		n2 = getRandomProtein();
		p1 = (Protein*) (*vertexes)[n1];
		p2 = (Protein*) (*vertexes)[n2];

		tries--;
	}


	if (tries <= 0) {

		TRACE(3,"Failed mutation. Reason: no uncomplexed proteins found.");
		return -1;
	}

	TRACE(3,"First Protein # " << p1->getMoleculeID() << " chosen.");
	TRACE(3,"Second Protein # " << p2->getMoleculeID() << " chosen.");

	/*
	 * Create new protein.
	 */
	Protein* newProtein = new Protein();

	newProtein->getConcentration()->value = DEFAULT_PROTEIN_CONC;

	// "P1" and "P2" (name+ID, not just name), for  example, would
	// complex to form the name "(P1_P2)"
	stringstream s;
	s << "(" << p1->name << p1->getMoleculeID() << "_" << p2->name << p2->getMoleculeID() << ")";
	newProtein->name = s.str();

	//get nodes for selected proteins
	ListDigraph::Node A = equations->nodeFromId(p1->nodeID);
	ListDigraph::Node B = equations->nodeFromId(p2->nodeID);

	//create node for new protein
	ListDigraph::Node AB = addProtein(newProtein);

	//create quantity variables
	Quantity* kDEG = new Quantity();
	Quantity* kFWD = new Quantity();
	Quantity* kREV = new Quantity();

	//set rates
	kDEG->value = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	kFWD->value = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	kREV->value = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);

	//add forward reactions
	ListDigraph::Arc A_AB = addInteraction(A, AB, kFWD, 1.0, 0.5,"fwd");
	ListDigraph::Arc B_AB = addInteraction(B, AB, kFWD, 1.0, 0.5,"fwd");

	//add reverse reactions
	ListDigraph::Arc AB_A = addInteraction(AB, A, kREV, 0.5, 1.0,"rev");
	ListDigraph::Arc AB_B = addInteraction(AB, B, kREV, 0.5, 1.0,"rev");

	ListDigraph::Arc AB_deg = addInteraction(AB, nullNode, kDEG, 1.0, 0,"deg");

	//TODO 6/4/11 document Term class
	Term* t1 = (*edges)[A_AB]->getTerm();
	t1->addFactor(p1->getConcentration());

	TRACE(3,"New Protein-Protein Complex in Cell " << (*edges)[A_AB]->name);

	Term* t2 = (*edges)[B_AB]->getTerm();
	t2->addFactor(p2->getConcentration());

	TRACE(3,"New Protein-Protein Complex in Cell " << cellID << " created.");

	return 1;
}

/**
 * Fill in
 */
int Cell::newBasic()
{

	TRACE(4,"Creating new Basic Protein in Cell " << cellID << "...");
	ListDigraph::Node A = addProtein(DEFAULT_PROTEIN_CONC);
	ListDigraph::Node B = addMRNA(DEFAULT_MRNA_CONC);
	ListDigraph::Node C = addDNA(DEFAULT_DNA_CONC);

	mutationLog << " New Basic Protein" << endl;;

	double kTXN = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	double kTSLN = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);

	double kDEGp = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	double kDEGmRNA = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);

	addInteraction(B, A, kTSLN, 0.0, 1.0, "tsln");
	addInteraction(C, B, kTXN, 0.0, 1.0, "txn");
	addInteraction(A,nullNode,kDEGp, 1.0, 0, "deg");
	addInteraction(B,nullNode,kDEGmRNA, 1.0, 0, "deg");

	TRACE(4,"Basic Protein in Cell " << cellID << " created.");
	return 1;
}

/**
 * Random rate change of a forward (protein -> ptm, protein -> complex) or translation/transcription interaction.
 */
int Cell::fwdRateChange()
{
	mutationLog << " Forward Rate Change" << endl;

	//temporary vector to hold forward rate interactions
	vector<Interaction*> fwdInteractions;

	//iterate through all interactions in the graph
	for (ListDigraph::ArcIt a(*equations); a != INVALID; ++a){

		//if this is a forward rate, or a transcription / translation rate, add it to the vector
		if( !strcmp((*edges)[a]->name, "fwd") || !strcmp((*edges)[a]->name, "txn") || !strcmp((*edges)[a]->name, "tsln")){
			fwdInteractions.push_back((*edges)[a]);

			TRACE(5, "adding " << (*vertexes)[equations->source(a)]->name << (*vertexes)[equations->source(a)]->getMoleculeID() << " -" << (*edges)[a]->getRate() << "-> " << (*vertexes)[equations->target(a)]->name << (*vertexes)[equations->target(a)]->getMoleculeID());
		}

	}

	//get a random number between 0 and the size of the forward interactions vector
	int num = r.randInt(fwdInteractions.size()-1);

	TRACE(5, "Forward interaction chosen for rate change");

	//get a new rate
	double newRate = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);

	//set the new rate for the chosen interaction
	fwdInteractions.at(num)->setRate(newRate);

	return 1;
}
/**
 * Random rate change of a reverse interaction.
 */
int Cell::revRateChange()
{
	mutationLog << " Reverse Rate Change" << endl;

	//hold forward rate interactions to randomly choose from
	vector<Interaction*> revInteractions;

	//iterate through all interactions in the graph
	for (ListDigraph::ArcIt a(*equations); a != INVALID; ++a){

		//if this is a reverse rate
		if( !strcmp((*edges)[a]->name, "rev")){
			//add it to the vector
			revInteractions.push_back((*edges)[a]);

			TRACE(5, "adding " << (*vertexes)[equations->source(a)]->name << (*vertexes)[equations->source(a)]->getMoleculeID() << " -" << (*edges)[a]->getRate() << "-> " << (*vertexes)[equations->target(a)]->name << (*vertexes)[equations->target(a)]->getMoleculeID());
		}

	}
	if(revInteractions.size() == 0)
	{
		mutationLog << " FAILED: no reverse rates exist\n";
		TRACE(5, "No reverse interactions exist");
		return 0;
	}

	//get a random number between 0 and the size of the forward interactions vector
	int num = r.randInt(revInteractions.size()-1);

	TRACE(5, "Reverse interaction chosen for rate change");

	//get a new rate
	double newRate = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);

	//set the new rate for the chosen interaction
	revInteractions.at(num)->setRate(newRate);


	return 1;
}

/**
 * Random rate change of a degradation interaction.
 */
int Cell::degRateChange()
{
	mutationLog << " Degradation Rate Change" << endl;

	//choose a random number between 0 and the number of degrading molecules
	int n = r.randInt(countInArcs(*equations, nullNode)-1);

	Interaction * selection;

	//iterate to the nth arc incoming to the null node (degrading)
	ListDigraph::InArcIt a (*equations, nullNode);
	while(n != 0){
		++a;
		n--;
	}
	selection = (*edges)[a];
	double oldRate = selection->getRate();

	//choose a new random rate
	double newRate = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);



	//set the selected arc to the new rate
	selection->setRate(newRate);

	TRACE(5, "Cell " << cellID << " Degradation Rate Change: " << (*vertexes)[equations->source(a)]->name << (*vertexes)[equations->source(a)]->getMoleculeID()
			<< " from " << oldRate << " to " << newRate);

	return 1;
}

/**
 * Apply a Post-Translational Modification to a randomly chosen Protein.
 */
int Cell::newPostTransMod()
{
	mutationLog << " New PTM" << endl;

	TRACE(4,"Adding new Post-Translational Modification (PTM) in Cell " << cellID << "...");

	/*
	 * Choose random protein that can receive a PTM. After (# of proteins / 3)
	 * selections of proteins that cannot receive additional PTMs,
	 * give up
	 */

	TRACE(5,"Choosing protein...");

	ListDigraph::Node n;
	Protein* p;

	//get a random protein
	n = getRandomProtein();
	p = (Protein*) (*vertexes)[n];


	TRACE(5,"Protein # " << p->getMoleculeID() << " chosen.");

	/*
	 * Create new protein.
	 */
	Protein* newProtein = new Protein();
	newProtein->getConcentration()->value = DEFAULT_PROTEIN_CONC;
	newProtein->name = "PTM";

	/*
	 * Copy post-translational modifications to new protein.
	 */
	int arrayLength = sizeof(p->postTransMods)/sizeof(*(p->postTransMods));
	for (int i = 0; i < arrayLength; i++)
		newProtein->postTransMods[i] = p->postTransMods[i];

	/*
	 * Increase the count of a randomly chosen PTM in new protein.
	 */
	int postTransModIndex = r.randInt(3);
	newProtein->postTransMods[postTransModIndex] += 1;

	//get node for selected protein
	ListDigraph::Node A = equations->nodeFromId(p->nodeID);

	//create node for new protein
	ListDigraph::Node B = addProtein(newProtein);

	//set rates
	double kFWD = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	double kREV = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
	double kDEG = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);

	//add forward and reverse reactions
	ListDigraph::Arc A_B = addInteraction(A, B, kFWD, 1, 1, "fwd");
	ListDigraph::Arc B_A = addInteraction(B, A, kREV, 1, 1, "rev");
	ListDigraph::Arc B_N = addInteraction(B, nullNode, kDEG, 1, 0, "deg");



	TRACE(4,"New PTM in Cell " << cellID << " added.");

	return 1;
}

//TODO 6/4/11 implement use of specific PTM's.
//TODO 6/4/11 explain why this is different from a Forward Rate Change
/**
 * Apply a histone modification to a random transcription interaction. 
 */
int Cell::histoneMod() {

	mutationLog << " Histone Modification" << endl;

	//hold transcription rate interactions to randomly choose from
	vector<Interaction*> txnInteractions;

	//iterate through all interactions in the graph
	for (ListDigraph::ArcIt a(*equations); a != INVALID; ++a){

		//if this is a transcription rate
		if( !strcmp((*edges)[a]->name, "txn")){
			//add it to the vector
			txnInteractions.push_back((*edges)[a]);

			TRACE(5, "adding " << (*vertexes)[equations->source(a)]->name << (*vertexes)[equations->source(a)]->getMoleculeID() << " -" << (*edges)[a]->getRate() << "-> " << (*vertexes)[equations->target(a)]->name << (*vertexes)[equations->target(a)]->getMoleculeID());
		}
	}
	//get a random number between 0 and the number of translation interactions
	int num = r.randInt(txnInteractions.size()-1);

	//get a factor [0,1) repression, (1,2] activation
	double txnFactor = r.rand(2.0);

	//set the product factor of the chosen transcription interaction
	txnInteractions.at(num)->setProductFactor(txnFactor);

	return 1;


}


/**
 * Not implemented
 */
int Cell::newPromoter() {
	//placeholder message
	mutationLog << " New Protein-promoter (not implemented)" << endl;
	
	ListDigraph::Node n1 = getRandomProtein();
	Protein* p1 = (Protein*) (*vertexes)[n1];

	vector<Molecule*> dna;
	for(ListDigraph::NodeIt iter(*equations); iter != INVALID; ++iter)
	  if( (*vertexes)[iter]->name.compare("DNA") == 0)
	     dna.push_back((*vertexes)[iter]);

	int num = r.randInt(dna.size()-1);

	DNA* selected = (DNA*) dna.at(num);	
	ListDigraph::Node d = equations->nodeFromId(dna.at(num)->nodeID);

	mutationLog << p1->getMoleculeID() << " " << p1->name << endl;
	mutationLog << selected->getMoleculeID() << " " << selected->name << endl;
	int found = 0;
	for(ListDigraph::InArcIt a(*equations, d); a != INVALID; ++a){
		if( equations->oppositeNode(d,a) == n1)
			found = 1;
	}

	if(!found){
		double kf = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);
		double kr = KINETIC_RATE_LOWER_LIMIT + r.rand(KINETIC_RATE_UPPER_LIMIT - KINETIC_RATE_LOWER_LIMIT);

		addInteraction(n1, d, kf-kr, 1, 0, "pro");
		selected->setHill(2);
		selected->setBindRates(kf,kr);
	}
	else
		mutationLog << "Promoter interaction already exists" << endl;
	return -1;
}
